<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>dbdsdc.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="DBDSDC.1"></a><a href="dbdsdc.f.html#DBDSDC.1">DBDSDC</a>( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ,
     $                   WORK, IWORK, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      CHARACTER          COMPQ, UPLO
      INTEGER            INFO, LDU, LDVT, N
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      INTEGER            IQ( * ), IWORK( * )
      DOUBLE PRECISION   D( * ), E( * ), Q( * ), U( LDU, * ),
     $                   VT( LDVT, * ), WORK( * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="DBDSDC.21"></a><a href="dbdsdc.f.html#DBDSDC.1">DBDSDC</a> computes the singular value decomposition (SVD) of a real
</span><span class="comment">*</span><span class="comment">  N-by-N (upper or lower) bidiagonal matrix B:  B = U * S * VT,
</span><span class="comment">*</span><span class="comment">  using a divide and conquer method, where S is a diagonal matrix
</span><span class="comment">*</span><span class="comment">  with non-negative diagonal elements (the singular values of B), and
</span><span class="comment">*</span><span class="comment">  U and VT are orthogonal matrices of left and right singular vectors,
</span><span class="comment">*</span><span class="comment">  respectively. <a name="DBDSDC.26"></a><a href="dbdsdc.f.html#DBDSDC.1">DBDSDC</a> can be used to compute all singular values,
</span><span class="comment">*</span><span class="comment">  and optionally, singular vectors or singular vectors in compact form.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  This code makes very mild assumptions about floating point
</span><span class="comment">*</span><span class="comment">  arithmetic. It will work on machines with a guard digit in
</span><span class="comment">*</span><span class="comment">  add/subtract, or on those binary machines without guard digits
</span><span class="comment">*</span><span class="comment">  which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
</span><span class="comment">*</span><span class="comment">  It could conceivably fail on hexadecimal or decimal machines
</span><span class="comment">*</span><span class="comment">  without guard digits, but we know of none.  See <a name="DLASD3.34"></a><a href="dlasd3.f.html#DLASD3.1">DLASD3</a> for details.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  The code currently calls <a name="DLASDQ.36"></a><a href="dlasdq.f.html#DLASDQ.1">DLASDQ</a> if singular values only are desired.
</span><span class="comment">*</span><span class="comment">  However, it can be slightly modified to compute singular values
</span><span class="comment">*</span><span class="comment">  using the divide and conquer method.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  UPLO    (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment">          = 'U':  B is upper bidiagonal.
</span><span class="comment">*</span><span class="comment">          = 'L':  B is lower bidiagonal.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  COMPQ   (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment">          Specifies whether singular vectors are to be computed
</span><span class="comment">*</span><span class="comment">          as follows:
</span><span class="comment">*</span><span class="comment">          = 'N':  Compute singular values only;
</span><span class="comment">*</span><span class="comment">          = 'P':  Compute singular values and compute singular
</span><span class="comment">*</span><span class="comment">                  vectors in compact form;
</span><span class="comment">*</span><span class="comment">          = 'I':  Compute singular values and singular vectors.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  N       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The order of the matrix B.  N &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  D       (input/output) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment">          On entry, the n diagonal elements of the bidiagonal matrix B.
</span><span class="comment">*</span><span class="comment">          On exit, if INFO=0, the singular values of B.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
</span><span class="comment">*</span><span class="comment">          On entry, the elements of E contain the offdiagonal
</span><span class="comment">*</span><span class="comment">          elements of the bidiagonal matrix whose SVD is desired.
</span><span class="comment">*</span><span class="comment">          On exit, E has been destroyed.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  U       (output) DOUBLE PRECISION array, dimension (LDU,N)
</span><span class="comment">*</span><span class="comment">          If  COMPQ = 'I', then:
</span><span class="comment">*</span><span class="comment">             On exit, if INFO = 0, U contains the left singular vectors
</span><span class="comment">*</span><span class="comment">             of the bidiagonal matrix.
</span><span class="comment">*</span><span class="comment">          For other values of COMPQ, U is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDU     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array U.  LDU &gt;= 1.
</span><span class="comment">*</span><span class="comment">          If singular vectors are desired, then LDU &gt;= max( 1, N ).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  VT      (output) DOUBLE PRECISION array, dimension (LDVT,N)
</span><span class="comment">*</span><span class="comment">          If  COMPQ = 'I', then:
</span><span class="comment">*</span><span class="comment">             On exit, if INFO = 0, VT' contains the right singular
</span><span class="comment">*</span><span class="comment">             vectors of the bidiagonal matrix.
</span><span class="comment">*</span><span class="comment">          For other values of COMPQ, VT is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDVT    (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array VT.  LDVT &gt;= 1.
</span><span class="comment">*</span><span class="comment">          If singular vectors are desired, then LDVT &gt;= max( 1, N ).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Q       (output) DOUBLE PRECISION array, dimension (LDQ)
</span><span class="comment">*</span><span class="comment">          If  COMPQ = 'P', then:
</span><span class="comment">*</span><span class="comment">             On exit, if INFO = 0, Q and IQ contain the left
</span><span class="comment">*</span><span class="comment">             and right singular vectors in a compact form,
</span><span class="comment">*</span><span class="comment">             requiring O(N log N) space instead of 2*N**2.
</span><span class="comment">*</span><span class="comment">             In particular, Q contains all the DOUBLE PRECISION data in
</span><span class="comment">*</span><span class="comment">             LDQ &gt;= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))
</span><span class="comment">*</span><span class="comment">             words of memory, where SMLSIZ is returned by <a name="ILAENV.94"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a> and
</span><span class="comment">*</span><span class="comment">             is equal to the maximum size of the subproblems at the
</span><span class="comment">*</span><span class="comment">             bottom of the computation tree (usually about 25).
</span><span class="comment">*</span><span class="comment">          For other values of COMPQ, Q is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IQ      (output) INTEGER array, dimension (LDIQ)
</span><span class="comment">*</span><span class="comment">          If  COMPQ = 'P', then:
</span><span class="comment">*</span><span class="comment">             On exit, if INFO = 0, Q and IQ contain the left
</span><span class="comment">*</span><span class="comment">             and right singular vectors in a compact form,
</span><span class="comment">*</span><span class="comment">             requiring O(N log N) space instead of 2*N**2.
</span><span class="comment">*</span><span class="comment">             In particular, IQ contains all INTEGER data in
</span><span class="comment">*</span><span class="comment">             LDIQ &gt;= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))
</span><span class="comment">*</span><span class="comment">             words of memory, where SMLSIZ is returned by <a name="ILAENV.106"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a> and
</span><span class="comment">*</span><span class="comment">             is equal to the maximum size of the subproblems at the
</span><span class="comment">*</span><span class="comment">             bottom of the computation tree (usually about 25).
</span><span class="comment">*</span><span class="comment">          For other values of COMPQ, IQ is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK    (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
</span><span class="comment">*</span><span class="comment">          If COMPQ = 'N' then LWORK &gt;= (4 * N).
</span><span class="comment">*</span><span class="comment">          If COMPQ = 'P' then LWORK &gt;= (6 * N).
</span><span class="comment">*</span><span class="comment">          If COMPQ = 'I' then LWORK &gt;= (3 * N**2 + 4 * N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IWORK   (workspace) INTEGER array, dimension (8*N)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0:  successful exit.
</span><span class="comment">*</span><span class="comment">          &lt; 0:  if INFO = -i, the i-th argument had an illegal value.
</span><span class="comment">*</span><span class="comment">          &gt; 0:  The algorithm failed to compute an singular value.
</span><span class="comment">*</span><span class="comment">                The update process of divide and conquer failed.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Further Details
</span><span class="comment">*</span><span class="comment">  ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Based on contributions by
</span><span class="comment">*</span><span class="comment">     Ming Gu and Huan Ren, Computer Science Division, University of
</span><span class="comment">*</span><span class="comment">     California at Berkeley, USA
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">  Changed dimension statement in comment describing E from (N) to
</span><span class="comment">*</span><span class="comment">  (N-1).  Sven, 17 Feb 05.
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      DOUBLE PRECISION   ZERO, ONE, TWO
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      INTEGER            DIFL, DIFR, GIVCOL, GIVNUM, GIVPTR, I, IC,
     $                   ICOMPQ, IERR, II, IS, IU, IUPLO, IVT, J, K, KK,
     $                   MLVL, NM1, NSIZE, PERM, POLES, QSTART, SMLSIZ,
     $                   SMLSZP, SQRE, START, WSTART, Z
      DOUBLE PRECISION   CS, EPS, ORGNRM, P, R, SN
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.148"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      INTEGER            <a name="ILAENV.149"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>
      DOUBLE PRECISION   <a name="DLAMCH.150"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANST.150"></a><a href="dlanst.f.html#DLANST.1">DLANST</a>
      EXTERNAL           <a name="LSAME.151"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="ILAENV.151"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>, <a name="DLAMCH.151"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANST.151"></a><a href="dlanst.f.html#DLANST.1">DLANST</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           DCOPY, <a name="DLARTG.154"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>, <a name="DLASCL.154"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>, <a name="DLASD0.154"></a><a href="dlasd0.f.html#DLASD0.1">DLASD0</a>, <a name="DLASDA.154"></a><a href="dlasda.f.html#DLASDA.1">DLASDA</a>, <a name="DLASDQ.154"></a><a href="dlasdq.f.html#DLASDQ.1">DLASDQ</a>,
     $                   <a name="DLASET.155"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>, <a name="DLASR.155"></a><a href="dlasr.f.html#DLASR.1">DLASR</a>, DSWAP, <a name="XERBLA.155"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, DBLE, INT, LOG, SIGN
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Test the input parameters.
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
<span class="comment">*</span><span class="comment">
</span>      IUPLO = 0
      IF( <a name="LSAME.167"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'U'</span> ) )
     $   IUPLO = 1
      IF( <a name="LSAME.169"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'L'</span> ) )
     $   IUPLO = 2
      IF( <a name="LSAME.171"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPQ, <span class="string">'N'</span> ) ) THEN
         ICOMPQ = 0
      ELSE IF( <a name="LSAME.173"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPQ, <span class="string">'P'</span> ) ) THEN
         ICOMPQ = 1
      ELSE IF( <a name="LSAME.175"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPQ, <span class="string">'I'</span> ) ) THEN
         ICOMPQ = 2
      ELSE
         ICOMPQ = -1
      END IF
      IF( IUPLO.EQ.0 ) THEN
         INFO = -1
      ELSE IF( ICOMPQ.LT.0 ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( ( LDU.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDU.LT.
     $         N ) ) ) THEN
         INFO = -7
      ELSE IF( ( LDVT.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDVT.LT.
     $         N ) ) ) THEN
         INFO = -9
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.194"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DBDSDC.194"></a><a href="dbdsdc.f.html#DBDSDC.1">DBDSDC</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.0 )
     $   RETURN
      SMLSIZ = <a name="ILAENV.202"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 9, <span class="string">'<a name="DBDSDC.202"></a><a href="dbdsdc.f.html#DBDSDC.1">DBDSDC</a>'</span>, <span class="string">' '</span>, 0, 0, 0, 0 )
      IF( N.EQ.1 ) THEN
         IF( ICOMPQ.EQ.1 ) THEN
            Q( 1 ) = SIGN( ONE, D( 1 ) )
            Q( 1+SMLSIZ*N ) = ONE
         ELSE IF( ICOMPQ.EQ.2 ) THEN
            U( 1, 1 ) = SIGN( ONE, D( 1 ) )
            VT( 1, 1 ) = ONE
         END IF
         D( 1 ) = ABS( D( 1 ) )
         RETURN
      END IF
      NM1 = N - 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If matrix lower bidiagonal, rotate to be upper bidiagonal
</span><span class="comment">*</span><span class="comment">     by applying Givens rotations on the left
</span><span class="comment">*</span><span class="comment">
</span>      WSTART = 1
      QSTART = 3
      IF( ICOMPQ.EQ.1 ) THEN
         CALL DCOPY( N, D, 1, Q( 1 ), 1 )
         CALL DCOPY( N-1, E, 1, Q( N+1 ), 1 )
      END IF
      IF( IUPLO.EQ.2 ) THEN
         QSTART = 5
         WSTART = 2*N - 1
         DO 10 I = 1, N - 1
            CALL <a name="DLARTG.229"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>( D( I ), E( I ), CS, SN, R )
            D( I ) = R
            E( I ) = SN*D( I+1 )
            D( I+1 ) = CS*D( I+1 )
            IF( ICOMPQ.EQ.1 ) THEN
               Q( I+2*N ) = CS
               Q( I+3*N ) = SN
            ELSE IF( ICOMPQ.EQ.2 ) THEN
               WORK( I ) = CS
               WORK( NM1+I ) = -SN
            END IF
   10    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If ICOMPQ = 0, use <a name="DLASDQ.243"></a><a href="dlasdq.f.html#DLASDQ.1">DLASDQ</a> to compute the singular values.
</span><span class="comment">*</span><span class="comment">
</span>      IF( ICOMPQ.EQ.0 ) THEN
         CALL <a name="DLASDQ.246"></a><a href="dlasdq.f.html#DLASDQ.1">DLASDQ</a>( <span class="string">'U'</span>, 0, N, 0, 0, 0, D, E, VT, LDVT, U, LDU, U,
     $                LDU, WORK( WSTART ), INFO )
         GO TO 40
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If N is smaller than the minimum divide size SMLSIZ, then solve
</span><span class="comment">*</span><span class="comment">     the problem with another solver.
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.LE.SMLSIZ ) THEN
         IF( ICOMPQ.EQ.2 ) THEN
            CALL <a name="DLASET.256"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'A'</span>, N, N, ZERO, ONE, U, LDU )
            CALL <a name="DLASET.257"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'A'</span>, N, N, ZERO, ONE, VT, LDVT )
            CALL <a name="DLASDQ.258"></a><a href="dlasdq.f.html#DLASDQ.1">DLASDQ</a>( <span class="string">'U'</span>, 0, N, N, N, 0, D, E, VT, LDVT, U, LDU, U,
     $                   LDU, WORK( WSTART ), INFO )
         ELSE IF( ICOMPQ.EQ.1 ) THEN
            IU = 1
            IVT = IU + N
            CALL <a name="DLASET.263"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'A'</span>, N, N, ZERO, ONE, Q( IU+( QSTART-1 )*N ),
     $                   N )
            CALL <a name="DLASET.265"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'A'</span>, N, N, ZERO, ONE, Q( IVT+( QSTART-1 )*N ),
     $                   N )
            CALL <a name="DLASDQ.267"></a><a href="dlasdq.f.html#DLASDQ.1">DLASDQ</a>( <span class="string">'U'</span>, 0, N, N, N, 0, D, E,
     $                   Q( IVT+( QSTART-1 )*N ), N,
     $                   Q( IU+( QSTART-1 )*N ), N,
     $                   Q( IU+( QSTART-1 )*N ), N, WORK( WSTART ),
     $                   INFO )
         END IF
         GO TO 40
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( ICOMPQ.EQ.2 ) THEN
         CALL <a name="DLASET.277"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'A'</span>, N, N, ZERO, ONE, U, LDU )
         CALL <a name="DLASET.278"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'A'</span>, N, N, ZERO, ONE, VT, LDVT )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Scale.
</span><span class="comment">*</span><span class="comment">
</span>      ORGNRM = <a name="DLANST.283"></a><a href="dlanst.f.html#DLANST.1">DLANST</a>( <span class="string">'M'</span>, N, D, E )
      IF( ORGNRM.EQ.ZERO )
     $   RETURN
      CALL <a name="DLASCL.286"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ORGNRM, ONE, N, 1, D, N, IERR )
      CALL <a name="DLASCL.287"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, IERR )
<span class="comment">*</span><span class="comment">
</span>      EPS = <a name="DLAMCH.289"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Epsilon'</span> )
<span class="comment">*</span><span class="comment">
</span>      MLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
      SMLSZP = SMLSIZ + 1
<span class="comment">*</span><span class="comment">
</span>      IF( ICOMPQ.EQ.1 ) THEN
         IU = 1
         IVT = 1 + SMLSIZ
         DIFL = IVT + SMLSZP
         DIFR = DIFL + MLVL
         Z = DIFR + MLVL*2
         IC = Z + MLVL
         IS = IC + 1
         POLES = IS + 1
         GIVNUM = POLES + 2*MLVL
<span class="comment">*</span><span class="comment">
</span>         K = 1
         GIVPTR = 2
         PERM = 3
         GIVCOL = PERM + MLVL
      END IF
<span class="comment">*</span><span class="comment">
</span>      DO 20 I = 1, N
         IF( ABS( D( I ) ).LT.EPS ) THEN
            D( I ) = SIGN( EPS, D( I ) )
         END IF
   20 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      START = 1
      SQRE = 0
<span class="comment">*</span><span class="comment">
</span>      DO 30 I = 1, NM1
         IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Subproblem found. First determine its size and then
</span><span class="comment">*</span><span class="comment">        apply divide and conquer on it.
</span><span class="comment">*</span><span class="comment">
</span>            IF( I.LT.NM1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        A subproblem with E(I) small for I &lt; NM1.
</span><span class="comment">*</span><span class="comment">
</span>               NSIZE = I - START + 1
            ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        A subproblem with E(NM1) not too small but I = NM1.
</span><span class="comment">*</span><span class="comment">
</span>               NSIZE = N - START + 1
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        A subproblem with E(NM1) small. This implies an
</span><span class="comment">*</span><span class="comment">        1-by-1 subproblem at D(N). Solve this 1-by-1 problem
</span><span class="comment">*</span><span class="comment">        first.
</span><span class="comment">*</span><span class="comment">
</span>               NSIZE = I - START + 1
               IF( ICOMPQ.EQ.2 ) THEN
                  U( N, N ) = SIGN( ONE, D( N ) )
                  VT( N, N ) = ONE
               ELSE IF( ICOMPQ.EQ.1 ) THEN
                  Q( N+( QSTART-1 )*N ) = SIGN( ONE, D( N ) )
                  Q( N+( SMLSIZ+QSTART-1 )*N ) = ONE
               END IF
               D( N ) = ABS( D( N ) )
            END IF
            IF( ICOMPQ.EQ.2 ) THEN
               CALL <a name="DLASD0.353"></a><a href="dlasd0.f.html#DLASD0.1">DLASD0</a>( NSIZE, SQRE, D( START ), E( START ),
     $                      U( START, START ), LDU, VT( START, START ),
     $                      LDVT, SMLSIZ, IWORK, WORK( WSTART ), INFO )
            ELSE
               CALL <a name="DLASDA.357"></a><a href="dlasda.f.html#DLASDA.1">DLASDA</a>( ICOMPQ, SMLSIZ, NSIZE, SQRE, D( START ),
     $                      E( START ), Q( START+( IU+QSTART-2 )*N ), N,
     $                      Q( START+( IVT+QSTART-2 )*N ),
     $                      IQ( START+K*N ), Q( START+( DIFL+QSTART-2 )*
     $                      N ), Q( START+( DIFR+QSTART-2 )*N ),
     $                      Q( START+( Z+QSTART-2 )*N ),
     $                      Q( START+( POLES+QSTART-2 )*N ),
     $                      IQ( START+GIVPTR*N ), IQ( START+GIVCOL*N ),
     $                      N, IQ( START+PERM*N ),
     $                      Q( START+( GIVNUM+QSTART-2 )*N ),
     $                      Q( START+( IC+QSTART-2 )*N ),
     $                      Q( START+( IS+QSTART-2 )*N ),
     $                      WORK( WSTART ), IWORK, INFO )
               IF( INFO.NE.0 ) THEN
                  RETURN
               END IF
            END IF
            START = I + 1
         END IF
   30 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Unscale
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="DLASCL.380"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ONE, ORGNRM, N, 1, D, N, IERR )
   40 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Use Selection Sort to minimize swaps of singular vectors
</span><span class="comment">*</span><span class="comment">
</span>      DO 60 II = 2, N
         I = II - 1
         KK = I
         P = D( I )
         DO 50 J = II, N
            IF( D( J ).GT.P ) THEN
               KK = J
               P = D( J )
            END IF
   50    CONTINUE
         IF( KK.NE.I ) THEN
            D( KK ) = D( I )
            D( I ) = P
            IF( ICOMPQ.EQ.1 ) THEN
               IQ( I ) = KK
            ELSE IF( ICOMPQ.EQ.2 ) THEN
               CALL DSWAP( N, U( 1, I ), 1, U( 1, KK ), 1 )
               CALL DSWAP( N, VT( I, 1 ), LDVT, VT( KK, 1 ), LDVT )
            END IF
         ELSE IF( ICOMPQ.EQ.1 ) THEN
            IQ( I ) = I
         END IF
   60 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO
</span><span class="comment">*</span><span class="comment">
</span>      IF( ICOMPQ.EQ.1 ) THEN
         IF( IUPLO.EQ.1 ) THEN
            IQ( N ) = 1
         ELSE
            IQ( N ) = 0
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If B is lower bidiagonal, update U by those Givens rotations
</span><span class="comment">*</span><span class="comment">     which rotated B to be upper bidiagonal
</span><span class="comment">*</span><span class="comment">
</span>      IF( ( IUPLO.EQ.2 ) .AND. ( ICOMPQ.EQ.2 ) )
     $   CALL <a name="DLASR.423"></a><a href="dlasr.f.html#DLASR.1">DLASR</a>( <span class="string">'L'</span>, <span class="string">'V'</span>, <span class="string">'B'</span>, N, N, WORK( 1 ), WORK( N ), U, LDU )
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="DBDSDC.427"></a><a href="dbdsdc.f.html#DBDSDC.1">DBDSDC</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
